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The mathematical structure of the widely popular Sudoku puzzles is akin to typical hard constraint 
satisfaction problems that lie at the heart of many applications, including protein folding and the 
general problem of finding the ground state of a glassy spin system. Via an exact mapping of 
Sudoku into a deterministic, continuous-time dynamical system, here we show that the difficulty of 
Sudoku translates into transient chaotic behavior exhibited by the dynamical system. In particular, 
we show that the escape rate At, an invariant characteristic of transient chaos, provides a single 
scalar measure of the puzzle's hardness, which correlates well with human difficulty level ratings. 
Accordingly, 77 = —Xog^QK can be used to define a "Richter"-type scale for puzzle hardness, with 
easy puzzles falling in the range < 77 < 1, medium ones within 1 < 77 < 2, hard in 2 < 77 < 3 and 
ultra-hard with 77 > 3. To our best knowledge, there are no known puzzles with 77 > 4. 



In Sudoku, considered as one of the world's most pop- 
ular puzzles [1 , we have to fill in the cells of a 9 x 9 grid 
with integers 1 to 9 such that in all rows, all columns 
and in nine 3x3 blocks every digit appears exactly once, 
while respecting a set of previously given digits in some 
of the cells (the so-called clues). Sudoku is an exact cover 
type constraint satisfaction problem [2] and it is one of 
Karp's 21 NP-complete problems [3^, when generalized to 
N X N grids [4 . NP-complete problems are "intractable" 
(unless P=NP) [2] |5] in the sense that all known algo- 
rithms that compute solutions to them do so in expo- 
nential worst-case time (in the number of variables N); 
in spite of the fact that if given a candidate solution, it 
takes only polynomial time to check its correctness. 

The intractability of NP-complete problems has im- 
portant consequences, ranging from public-key cryptog- 
raphy to statistical mechanics. In the latter case, for the 
ground-state problem of Ising spin glasses (±1 spins), one 
needs to find the lowest energy configuration among all 
the 2^ possible spin configurations. Additionally, to de- 
scribe the statistical behavior of such Ising spin models, 
one has to compute the partition function, which is a sum 
over all the 2^ configurations. Barahona [6 , then Istrail 
[7] have shown that for non-planar crystalline lattices, the 
ground-state problem and computing the partition func- 
tion are NP-complete [7 . Since there is little hope in 
providing polynomial time algorithms for NP-complete 
problems, the focus shifted towards understanding the 
nature of the complexity forbidding fast solutions to these 
problems. There has been considerable work in this di- 
rection, especially for the Boolean satisfiability problem 
/c-SAT, which is NP-complete for k > 3. Due to com- 
pleteness, all problems in NP (hence Sudoku as well), 
can be translated (in polynomial time) and formulated 
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as a /c-SAT problem. In /c-SAT we are given N Boolean 
variables to which we need to assign Os or Is (TRUE or 
FALSE) such that a given set of clauses in conjunctive 
normal form are all satisfied (evaluate to TRUE). Just as 
for the spin glass model, here we also have exponentially 
many (2^) configurations or assignments to search. 

In the following we treat algorithms as dynamical sys- 
tems. An algorithm is a finite set of instructions acting in 
some state space, applied iteratively from an initial state 
until an end state is reached. For example, the simplest 
algorithm for the Ising model ground state problem, or 
the 3- SAT problem would be exhaustively testing poten- 
tially all the 2^ configurations, which quickly becomes 
forbidding with increasing N. To improve performance, 
algorithms have become more sophisticated by exploit- 
ing the structure of the problem (of the state space). Ac- 
cordingly, now 3- SAT can be solved by a deterministic 
algorithm with an upper bound of 0(1.473^) steps [8]. 
Here we will only deal with deterministic algorithms that 
is, once an initial state is given, the "trajectory" of the 
dynamical system is uniquely determined. Thus, we ex- 
pect that the dynamics of those algorithms that exploit 
the structure of hard problems will reflect the complex- 
ity inherent in the problem itself. Complex behavior by 
deterministic dynamical systems is coined chaos in the 
literature [9HTT] , and thus the behavior of algorithms for 
hard problems is expected to appear highly irregular or 
chaotic [12 . 

Although the theory of nonlinear dynamical systems 
and chaos is well-established, it has not yet been ex- 
ploited in the context of optimization algorithms. One of 
the difficulties lies with the fact that most optimization 
algorithms are discrete and not easily cast in forms 
amenable to chaos theory methods. Recently, however, 
we have provided [13 a deterministic continuous-time 
solver for the Boolean satisfiability problem /c-SAT using 
coupled ordinary differential equations (ODE) with a 
one-to-one correspondence between the /c-SAT solution 
clusters and the attractors of the corresponding system of 
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FIG. 1: Sudoku and its Boolean representation, (a) a typical puzzle with bold digits as clues (givens). (b) Setup of the 
Boolean representation in a 9 x 9 x 9 grid, (c) Layer L4 of the puzzle (the one containing the digit 4) with 1-s in the location 
of the clues and the regions blocked out for digit 4 by the presence of the clues (shaded area). 



ODEs. This continuous-time dynamical system (CTDS) 
is in a form naturally suited for chaos theory methods, 
and thus it allows us to study the relationship between 
optimization hardness and chaotic behavior. Here we 
will focus only on solvable (SATisfiable) instances, and 
thus the observed chaotic behavior will necessarily be 
transient [TTJ [Ml [15] . We need to emphasize, however, 
that the dynamical properties characterize both the 
problem and the algorithm itself. For this reason, one 
compares the dynamical properties across problems of 
varying hardness using the same algorithm. Neverthe- 
less, since there are problem instances that are hard for 
all known algorithms, the appearance of transient chaos 
should be a universal feature of hard problems. It is 
also important to observe that transient chaos is not an 
N ^ 00 asymptotic behavior, but it appears for finite 
and thus measures of chaos can be used to characterize 
and categorize the hardness of individual instances of 
finite problems. To illustrate this, here we first map 
the popular 9x9 (hence finite) version of Sudoku 
into /c-SAT, then we solve it using our deterministic 
continuous-time solver [13]. By analyzing the behavior 
of the corresponding trajectories of the CTDS we show 
the appearance of transient chaos when increasing the 
hardness of the Sudoku problems, and show that the 
level of hardness (taken from human ratings of the puz- 
zles) correlates well with a chaotic invariant, namely the 
lifetime of chaos where is called the escape rate 
fll . We conclude with a discussion on algorithmic per- 
formance, dynamical properties and problem complexity. 



Results 

Sudoku as k-SAT 

Because our continuous-time dynamical system [13] 
was designed to solve /c-SAT formulae in conjunctive nor- 



mal form (CNF), we first briefly describe how Sudoku can 
be interpreted as a +l-in-9-SAT formula, and then how 
it is transformed into the standard CNF form. Further 
details are shown in the Methods section. 

In a Sudoku puzzle we are given a square grid with 
9x9 = 81 cells, each to be filled with one of nine symbols 
(digits) Dij G {1, . . . , 9}, z, j = 1, . . . , 9 (with the upper- 
left corner of the puzzle corresponding to i = 1, j = 1). 
When the puzzle is completed each of the columns, rows 
and 3x3 sub-grids (blocks partitioned by bold lines. Fig. 
[1^) must contain all the 9 symbols. Equivalently, all 9 
symbols must appear once and only once in each row, 
column and 3x3 sub-grid. 

To formulate Sudoku as a constraint satisfaction prob- 
lem (CSP) using Boolean variables, we associate to each 
symbol (digit) an ordered set of 9 Boolean variables 
(TRUE="1", FALSE="0"). The digit Aj in cell (ij) 
will be represented as the ordered set {xjj^ . . . , x^^) with 
x^^- G {0, 1}, a = 1, . . . , 9, such that always one and only 
one of them is 1 (TRUE). Thus Dij = a is equivalent to 
writing x^j = ^^,6, where Sa,b is the Kronecker delta func- 
tion. This way we have in total 9 x 9 x 9 = 729 Boolean 
variables xfj^ which we can picture as being placed on a 
3D grid (Fig. [TJ)), with a corresponding to the grid in- 
dex along the vertical direction, and hence a is the digit 
that is filling the corresponding (i, j) cell in the original 
puzzle. The corresponding 9x9 2D layer at height a will 
be denoted by L^. For example in the puzzle shown in 
Fig. ^ Di^g = 4. In the given vertical column the vari- 
able in the a*^ cell is xf g = (5^,4- The Sudoku constraints 
can also be simply encoded using Boolean variables (see 
Methods). They come from: 1) uniqueness of the sym- 
bols in all the (i, j) Sudoku cells, 2) a symbol must occur 
once and only once in each row, column and in each of 
the nine 3x3 subgrids, and 3) obeying the clues. Con- 
straint type 1) was already expressed above, namely that 
for every cell (i, j), in the set {xjj^ . . . , x^^) one and only 
one variable is TRUE, all others must be FALSE. Type 
2) constraints are similar, e.g., in row i and layer a the 
set {xfi^ . . . ,^^9) must contain one and only one TRUE 
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variable, all others must be false and this must hold for 
all rows and layers, etc. Observe that all constraints are 
in the form of a set of 9 Boolean variables of which we 
demand that one and only one of them be TRUE, all 
others FALSE. When this is satisfied, we say that the 
constraint itself (or "clause" ) is satisfied, or TRUE. Such 
CSPs are called +l-in-/c-SAT and they are part of so- 
called "locked occupation problems", which is a class of 
exceptionally hard CSPs [16l ^7\. Type 3) constraints 
are generated by the clues (or givens) which are symbols 
already filled in some of the cells and their number and 
positioning determines the difficulty of the puzzle. They 
are also set in a way to guarantee a unique solution to 
the whole puzzle. If there are given d clues, then this im- 
plies setting d Boolean variables to TRUE, which means 
eliminating exactly M constraints of type 1) and 2) (one 
vertical or uniqueness constraint, one row, one column 
and one 3x3 subgrid constraint). Thus, Sudoku is a 
+l-in-9-SAT type CSP with N Boolean variables and 
324 — Ad constraints. N is di complicated function of the 
positioning of the clues. 

In order to apply our continuous-time SAT solver we 
need to bring the +l-in-9-SAT type CSP above into con- 
junctive normal form. In /c-SAT there are N Boolean 
variables Xi = {0, 1} and an instance is given as a 
propositional formula J^, which is the conjunction (AND, 
denoted by A) of M clauses (constraints) Cm- ^ = 
Ci A • • • A Cm A • • • A Cm • Each clause is the disjunction 
(OR, denoted by V) of k literals. A literal is a variable 
(xi) or its negation {xi). For example a 3- SAT constraint 
could be Ci = xi V x^V x^. All Boolean propositions T 
can be formulated in CNF. 

Once the transformation to CNF is completed we are 
left with N variables and M SAT clauses (see Methods). 
We will denote the number of variables appearing in con- 
straint m by km-, = 1, . . . , M (clearly, 1 < /Cm ^ 9). 
The parameters M and {km}m=i depend on the 
clues that are difficult to express analytically, but easy to 
determine computationally, as illustrated via examples. 



The continuous-time deterministic A;-SAT solver 

In Ref [1^ a continuous-time deterministic solver was 
introduced to solve /c-SAT problems in conjunctive nor- 
mal form. The set of clauses specifying the constraints 
are translated into an M x matrix: C = {cmi} with 
Cmi = 1 if the variable Xi is present in clause m in di- 
rect (non- negated) form, namely Xi G C^, Cmi = — 1 if 
Xi G Cm and Cmi = if and Xi are both absent from 
Cm- To every variable Xi one associates a continuous 
spin variable Si G [—1, 1] such that when = ±1 then 
Xi = (1 ^ Si)/2 G {0,1}, and to every clause Cm one 
associates the function: 

N 

Km{s) = 2-^-\{{l-CmjSj) , mG {!,..., M}. (1) 



We have Km e [0, 1] for ah s G [-1, 1]^. It is easy to 
check that Km = only for those Si G { — 1,+!} values 
for which the corresponding Xi-s satisfy clause Cm (oth- 
erwise we always have Km > 0) . That is. Km plays the 
role of an energy function for clause Cm and its ground 
state value of Km = is reached if Cm is TRUE, and only 
then. We also need the quantities Kmi = Km/ {1 — CmiSi) 
that is, with the i-th term missing from the product in 
([T]). Clearly, Kmi ^ [O^l/^]. The continuous time dy- 
namical system introduced in [13 is defined via the set 
of (A^ -h M) ordinary differential equations (ODEs): 



dsi 
dt 

dOjm 

~dt 



M 



^2amCmiKmi{s)Km{s), i = l,...,N (2) 



amKm{s), m = 1, . . . ,M , 



(3) 



with the only requirements that 5^(0) G [—1,1], Vi and 
<^m(0) > 0, Vm. The latter implies from ([3| that 
<^m(^) > 0, Vm, t. It was shown in Ref [13 that system 
([2][3| always finds the solutions to /c-SAT problems (en- 
coded via the C matrix), when they exist, from almost all 
initial conditions (the exception being a set of Lebesgue 
measure zero). Here we give an intuitive picture for why 
that is the case. Due to (|3| the auxiliary variables am 
grow exponentially at rate Km- That is, the further is 
Km from its ground-state value of 0, the faster am grows 
(in that instant). Moreover, the longer has Km been 
away from zero, the larger is a^, as seen from the formal 

solution to (3): am{t) = am{0) exp (^J^ drKm^ - Equa- 
tion ([2| can equivalently be written as a gradient descent 
on an energy landscape V(s,a), that is ds/dt = — VgF, 
where Vs is the gradient operator in the spin variables 
and V{s,a) = X]m^^^m(^) • Clearly, V > \/t and 
V = if and only if s is a /c-SAT solution, i.e., satisfies 
all the clauses {Km{s) = 0, Vm). From the behavior of 
the am variables discussed above it also follows that the 
least satisfied constraints will dominate V (terms with 
the largest a^-s). Without restricting generality, let the 
aiKf term be the most dominant at t. Then keeping only 
the dominant term on the rhs of ([2| for those i for which 
cii 7^ we get dsi/dt = 2aiCii{l — ciiSi){Ki 
alently: d{l — ciiSi)/dt = —(1 — CiiSi)2ai{Kii) 
shows that the term {l—cuSi) is driven exponentially fast 
towards zero, that is towards satisfying Ki (and all the 
other constraints containing this term). As Ki decreases, 
some other constraint becomes dominant, and thus, in 
a continuous fashion, all constraints are driven towards 
satisfiability. The exponential growth guarantees that 
the trajectory is always pulled out of any potential well. 
When the problem is unsatisfiable, the system generates 
a chaotic dynamics in [—1,1]^, indefinitely. For more 
details about the properties of the CTDS ([2][3| see Ref. 
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FIG. 2: Solving Sudoku puzzles with the deterministic continuous- time solver (2] 3). (a) presents an easy puzzle 
with the evolution of the continuous-time dynamics shown within a 3 x 3 grid (rows 4-6, columns 7-9). (b) shows the same, 
but for a known, very hard puzzle called Platinum Blonde ^19^ . 



Puzzle hardness as transient chaotic dynamics 

Since Sudoku puzzles always have a solution, the cor- 
responding Boolean SAT CNF formulation also has a so- 
lution, and system (2]|3) will always find it. The nature 
of the dynamics, however will depend on the hardness of 
the puzzle as we describe next. 

In Figj2^ we show an easy puzzle with 34 clues (black 
numbers) [22 . After transforming this problem into 
SAT, we obtain = 126 and M = 717, with a constraint 
density of a = M/N = 5.69. As described above, in our 
implementation there is a spin variable s'^j associated to 
every Boolean variable oc^ in every 3D cell {i^j^a). In 
the right panels of Fig. [2] we show the dynamics of the 
spin variables in the cells of the 3x3 grid formed by 
rows 4-6 and columns 7-9. The s^j{t) curves are colored 
by the digit a they represent (a = 1, . . . , 9) as indicated 
in the color legend of Fig |2j The dynamics was started 
from a random initial condition. Indeed, our solver finds 
the solution very quickly, for the easy puzzle in Fig|2^. 

In Fig. we show the dynamical evolution of vari- 
ables for a very hard Sudoku instance with only 21 clues. 
This puzzle has been listed as one of the world's hard- 



est Sudokus, and even has a special name: "Platinum 
Blonde" [H [19], and it was the most "difficult" for our 
solver among all the puzzles we tried. After transform- 
ing it into SAT CNF, we obtain = 257 variables and 
M = 2085 constraints. Not only that we have twice 
as many unknown variables but the constraint density 
a = M/N — 8.11 is also larger than in the previous 
case, signaling the hardness of the corresponding SAT 
instance. The complexity of the dynamics in this case is 
seen in the right panel of Fig. [2]3, exhibiting long chaotic 
transients before the solution is found at around t ^ 150. 
For an animation of the dynamics for a similarly hard 
puzzle [12 see Ref [20j. 

We can also observe from the right panels in Fig |2] that 
there is one dominating digit (a- value), corresponding to 
which vertical cell at that given (i, j) grid cell has the 
largest \sfj \ value. This can be taken as the digit Dij the 
solver is considering in the given grid cell (i, j) at that 
moment. We will use this observation to provide below an 
alternate illustration of the dynamics' transiently chaotic 
behavior. Let us fix a random initial condition except for 
two chosen variables that are varied along the points of 
a square grid within the domain [—1,1]^. There is no 
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t = 10 t = 15 t = 20 




S2 S2 S2 

FIG. 3: Puzzle hardness as chaotic dynamics. We color the points of a 10*^ x 10^ grid in an arbitrary plane (si,S2) at 
time instant t according to the digit Dpq the solver is considering in an arbitrary but fixed cell (p, q) at that instant, given that 
we started the trajectory of the CTDS from those grid-points. For these initial conditions only the points in the (si, S2) plane 
were varied, all other spin values were kept fixed at the same randomly chosen values. For an easy problem (top row of panels), 
and for {p^q) = (1, 1) almost all initial conditions in this plane involve only two digits, and after t = 20 the corresponding 
trajectories have converged to the solution digit (9, light blue), except for a thin line, which, however, will also become light 
blue. The bottom row of panels shows the same for a hard problem based on what happens in the cell {p^q) — (6,8). The 
strong sensitivity to initial conditions appears as fractal structures of increasing complexity as time goes on, before eventually 
everything converges to the same color/digit (not shown). 



particular relevance as to which pairs of variables are 
chosen to be varied, let us denote them by Si and 52. 
Let us choose an arbitrary empty cell (p, q) in the original 
Sudoku puzzle and monitor the dominating digit in it at 
time t. We will color the initial conditions in the plane 
(si, 52) according to the dominating digit in (p, q) at time 
t. This will provide a map expressing the "sensitivity 
to initial conditions" that varies across time. Since all 
puzzles have solutions, the maps eventually assume one 
solid color according to the digit of the solution in the 
monitored cell, however, for hard puzzles, it may assume 
highly complex patterns before it does that, as shown in 
Fig. [3] In Figj3]we show these colormaps for the easy 
and hard Sudoku puzzles shown in Figj2] at times t = 
10, 15, 20. For the easy puzzle (top row of panels) the cell 
was chosen to be (p, g) = (1,1). At time t = 10 the whole 
map shows Di i = 6 (orange), which is not the solution 
digit (it is still searching for the solution). At time t = 15, 
however, we see two clearly separated domains, in one of 



them Dii = 6, in the other Dn = 9 (cyan) and the latter 
is the correct digit. As time passes, the orange (incorrect) 
domain shrinks, because trajectories from an increasing 
number of initial conditions find the solution. At t = 20 
almost the whole map shows the correct digit Dn = 9, 
except for a thin line. 

In the case of the hard Sudoku puzzle (bottom row in 
Fig. [sj (p, q) = (6, 8)) more colors enter the picture with 
time, in a complex fractal-like pattern. On this fractal 
set changing the initial condition slightly may result in a 
completely different digit (color) being considered in cell 
(p, q) at time t. This sensitivity to initial conditions is 
indicative of the chaotic behavior of the (deterministic) 
search dynamics. 

The appearance of transient chaos is a fundamental 
feature of the search dynamics and can be used to sep- 
arate problems by their hardness. In Ref [13] we have 
shown that within the thermodynamic limit {N 00, 
M ^ 00, a = M/N — const.) of random /c-SAT ensem- 
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bles this appears as a phase transition at the so-cahed 
chaotic transition point in terms of the constraint den- 
sity a = M/N. Since there is no "thermodynamic hmit" 
for 9 X 9 Sudoku problems {N < 729), one cannot define a 
simple order-parameter and use it to rate problem hard- 
ness in the same way [13]. However, once a problem is 
given, the corresponding dynamical system ([2][3| is well 
defined, and so is its dynamical behavior. Even though 
we do not have a well-defined ensemble-based statistical 
order parameter, (which has little meaning for specific 
SAT instances anyway), here we show next how can we 
use a well-known invariant quantity from non-linear dy- 
namical system's theory to categorize problem hardness 
for specific instances. 

A Richter-type scale for Sudoku hardness 

As suggested by the two examples in Fig [3] the hard- 
ness of Sudoku puzzles correlates with the length of 
chaotic transients. A consistent way to characterize these 
chaotic transients is to plot the distribution of their life- 
time. Starting trajectories from many random initial 
conditions, let p{t) indicate the probability that the dy- 
namics has not found the solution by analog time t. A 
characteristic property of transient chaos JTJ |2l] in hy- 
perbolic dynamical systems is that p{t) shows an expo- 
nential decay: p{t) ~ e~'^^, where is called the escape 
rate. The escape rate, an easily measurable quantity, 
theoretically can be expressed as a zero of the spectral 
determinant of the evolution operator corresponding to 
the dynamical system ( [2p| ) and well approximated using 
the machinery of cycle expansions based on dynamical 
zeta functions [2j. It is an invariant measure of the dy- 
namics in the sense that it characterizes solely the chaotic 
non-attracting set in the phase space of the system, and 
it does not depend on the distribution of the initial condi- 
tions, its support, or the details of the region from where 
the escape is measured (as long as it contains the non- 
attracting set) [TT. 

In Fig. [4^ we plot the distribution p{t) in log-linear 
scale for several puzzles gathered form the literature. The 
distributions were obtained from over 10"^ random initial 
conditions. The decay shows a wide range of variation 
between the puzzles. For easy puzzles the transients are 
very short, p{t) decays fast resulting in large escape rates 
but for hard puzzles k can be very small. Fig. shows a 
zoom onto the p{t) of hard puzzles. In spite of the large 
variability of the decay rates, we see that in all cases 
the escape is exponentially fast or faster (the curves in 
Figures |4^,b are straight lines or bend downward). 

The several orders of magnitude variability of natu- 
rally behooves us to use a logarithmic measure of hi for 
puzzle hardness, see Figj4]3, which shows the escape rates 
on a semilog scale as function of the number of clues, d. 
Thus, the escape rate can be used to define a kind of 
"Richter"-type scale for Sudoku hardness: 

^ = -logio(^) (4) 



with easy puzzles falling in the range < 7^ < 1, medium 
ones in 1 < ?7 < 2, hard ones in 2 < 77 < 3 and for 
ultra-hard puzzles 77 > 3. 

We chose several instances from the "Sudoku of the 
Day" website [22 in four of the categories defined there: 
easy (black square), medium (red circle), hard (green 
x) and absurd (blue star). These ratings on the web- 
site try to estimate the hardness of puzzles when solved 
by humans. These ratings correlate very well with our 
hardness measure rj, giving an average hardness value of 
{t]) = 0.816 for easy, {r]) = 1.439 medium, {rj) = 1.782 for 
hard and {r]) = 1.809 for what they call absurd. Another 
site we analyzed puzzles from is "Extreme Sudoku" [23] 
(brown + signs on FigQ. It claims to offer extremely 
hard Sudoku puzzles, their categories being: evil, exces- 
sive, egregious, excruciating and extreme. Indeed those 
puzzles are difficult with a range of r] G [1.1, 1.9] on the 
hardness scale, however, still far from the hardest puzzles 
we have found in the literature. Occasionally, daily news- 
papers present puzzles claimed to be the hardest Sudoku 
puzzles of the year. In particular, the escape rate for the 
Caveman Circus 2009 winner [24 (turquoise diamond) 
and the Guardian 2010 hardest puzzle [25] (maroon di- 
amond) are indeed one order of magnitude smaller than 
the hardest puzzles on the daily Sudoku websites, placing 
them at 7/ = 2.93 and r] — 2.82 on the hardness scale. The 
USA Today 2006 hardest puzzle [26j, however, does not 
seem to be that hard for our algorithm having t\ = 2.17 
(magenta diamond). Eppstein [27 gives two Sudoku ex- 
amples (orange left-pointing triangles) while describing 
his algorithm, one with r] = 1.288 and a much harder 
one with rj = 2.017. Elser et at. [T2 present an ex- 
tremely hard Sudoku (black filled circle), which has an 
escape rate oi n = 0.0023 resulting in 7/ = 2.639. 

The smallest escape rates we have found are for the 
Sudokus listed as the hardest on Wikipedia [TS] [28] (red 
triangles). The five puzzles, which we tested are called 
Platinum Blonde, Golden Nugget, Red Dwarf, coly013 
and tarx0134. They have a hardness in the range 3 < 
T] < 3.6, the Platinum Blonde (shown in Figj2]3) being 
the hardest with rj = 3.5789 (corresponding to an escape 
rate of = 0.00026). 

While the escape rate correlates surprisingly well with 
human ratings of Sudoku hardness, it is natural to ex- 
pect a correlation with the number of clues, d. Indeed, 
as a general rule of thumb, the fewer clues are given, 
the harder the puzzle, however, this is not universally 
true [1]. Here we tested a few instances with minimal 
[29 , that is 17 clues and almost minimal 18 clues (or- 
ange filled circles) [30-32 . As seen from Figj4]3, these 
are actually easier (1.2 < rj < 2.4) than the hardest in- 
stances with more d = 21,22 clues. In Figj4]i we then 
plot the escape rate as function of the constraint density 
a = M/N ^ leading to practically the same conclusion. 
This is because the constraint density a is essentially lin- 
early correlated with the number of givens d^ as shown 
in Figj4^. The apparent non-monotonic behavior of puz- 
zle hardness with the number of givens, (or constraint 
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FIG. 4: Escape rate as ha rdness indicator, (a) sliows tlie distribution in log-linear scale of the fraction p(t) of 10^ randomly 
started trajectories of (2|3) that have not yet found a solution by analog time t for a number of Sudoku puzzles taken from the 
literature (see legend and text) with a wide range of human difficulty ratings. The escape rate is obtained from the best fit to 
the tail of the distributions, (b) is a magnification of (a) for hard puzzles, (c) and (d) show the escape rate n in semilog scale 
vs the number of clues d and constraint density a indicating good correlations with human ratings (color bands), (e) shows 
the relationship between the number of clues d and a for the puzzles considered. 



density) is due to the fact that hardness cannot simply 
be characterized by a global, static variable such as d or 
a, but it also depends on the positioning pattern of the 
clues, as also shown by concrete examples in Ref [Ij. 



Discussion 

Using the world of Sudoku puzzles, here we have pre- 
sented further evidence that optimization hardness trans- 
lates into complex dynamical behavior by an algorithm 
searching for solutions in an optimal fashion. Namely, 
there seems to be a trade-off between algorithmic per- 
formance and the complexity of the algorithm and/or its 



behavior. Simple, sequential search algorithms have a 
trivial description and simple dynamics, but an abysmal 
worst-case performance (2^), whereas algorithms that 
are among the best performers are complex in their de- 
scription (instruction-list) and/or behavior (dynamics). 
This happens because in order to improve performance, 
algorithms have to exploit the structure of the problem 
one way or another. As hard problems have complex 
structures, the dynamics of the algorithms should be in- 
dicative of the problem's hardness. However, as a word of 
caution, observing complex dynamics performed by some 
black-box algorithm does not necessarily imply problem 
hardness. For example, one could consider any arbitrary, 
but ergodic dynamical system with complex behavior in 
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the same state space as the problem's. Ergodicity guar- 
antees the algorithm to eventually visit all of the 2^ 
states, and hence to always find solutions. But its in- 
struction list would have no relevance to the problem 
itself (apart from the checking instructions to see if the 
new state satisfies the problem) and thus, it could take 
long times to find solutions even for problems that are 
otherwise easily solved by other algorithms. Hence, dy- 
namical properties can only be regarded as descriptors 
of problem hardness if they are generated by algorithms 
that: 1) exploit the structure of the state space of the 
problem and 2) they show similar or better performance 
compared to other algorithms on the same problems. 

The continuous-time dynamical system [13] ([2][3| as a 
deterministic algorithm does have these features: 1) the 
search happens on an energy landscape V = ^m^m 
that incorporates simultaneously all the constraints 
(problem structure) 2) it solves easy problems efficiently 
(polynomial time, both analog and discrete) and 3) it 
guarantees to find solutions to hard problems even for 
solvable cases where many other algorithms fail. Al- 
though it is not a polynomial cost algorithm, it seems 
to find solutions in continuous-time t that scales poly- 
nomially with N [T3|. These features and the fact that 
the algorithm is formulated as a deterministic dynami- 
cal system with continuous variables, allows us to apply 
the theory of nonlinear dynamical systems on CTDS ^ 
[3| to characterize the hardness of Boolean satisfiability 
problems. In particular, via the measurable escape rate 

or its negative log- value 7^, we can provide a single- 
scalar measure of hardness, well defined for any finite 
instance. We have illustrated this here on Sudoku puz- 
zles, but the analysis can be repeated on any other en- 
semble from NP. Having a mathematically well-defined 
number to characterize optimization hardness for spe- 
cific problems in NP provides more information than the 
polynomial/exponential-time solvability classification, or 
knowing what the constraint density a = M/N is (the 
latter being a non-dynamic/static measure). Moreover, 
within the framework of CTDS ([2]|3|, dynamical systems 
and chaos theory methods can now be brought forth to 
help develop a novel understanding of optimization hard- 
ness. 



Methods 

Here we continue to describe in detail how a Sudoku 
puzzle is transformed into a SAT problem in CNF. 

Type 1) constraints (main text) impose the uniqueness 
of the symbol Dij in a given cell, expressed as a +l-in- 
9-SAT constraint: 

{^ij 7 -^ij 7 • • • 7 ^ij ) • (5) 

Having 9x9 cells in the puzzle, this gives in total 81, 



+l-in-9-SAT constraints. 

Type 2) constraints on rows, columns and sub-grids 
further impose that in every layer La we have the follow- 
ing 27, +l-in-9-SAT constraints: 



Rows: 

^i25 • • • 7 ^ig)' = 1, . . . , 9 

Columns: 

{Xij , , • • • , Xgj) j = 1, . . . , 9 

Subgrids: 



(6) 
(7) 



/ ™(2 (X CL 

V'^m+l,n+l' '^m+l,n+25 •^m+l,n+35 

^m+2,n+l5 ^m+2,n+25 ^m+2,n+35 (§) 
^m+3,n+l5 ^m+3,n+25 ^m+3,n+3) • = 0, 3, 6 

Together with the 81 constraints of type 1) we thus have 
in total 9 X 27 + 81 = 324 constraints in +l-in-9-SAT 
form. 

Finally, type 3) constraints are imposed via d given 
digits or clues. It was only recently shown that unique- 
ness of a solution demands that d > 17 [29]. As discussed 
in the main text, each clue will eliminate 4 constraints: 
in its vertical tower, its column, its row and the 3x3 
sub-grid containing the clue. For example, let us exam- 
ine layer L4 (Figjl]^) of the puzzle shown in Figjl^. There 
are three clues of 4 in cells (1,9), (3,3), (4,4) and thus 



^1,9 



1, X' 



- 1 ^4 _ 



5 -^3,3 — ^4,4 



1 have to be fixed as TRUE in 



L4. In order to satisfy the constraints, the other variables 
in the same rows, columns, blocks and vertical columns 
must be set to FALSE. The unknown variables left in 
the SAT problem will be those in the light cells of FigjT]^. 
(The other clues will eliminate constraints and variables 
in other layers and vertical columns.) The total num- 
ber of unknown variables N depends on d and on the 
placement of clues. The number of constraints is always 
324 — 4(i, however the number of variables in a clause can 
vary. For example in FigjT]^ the constraint correspond- 
ing to the second row in layer L4 has only 2 unknown 
variables left (+l-in-2-SAT). 

After the unknown Boolean variables and the con- 
straints have been identified we need to transform the 
formula into CNF. There are several ways of doing this, 
here we use the following general procedure. A -\-l-m-k- 
SAT clause defined on the (^i, ^2, • • • , ^fe) variables can 
be written as one A:-SAT and k{k - l)/2 of 2-SAT con- 
straints: 



(^1 V^2 V--- V^fc) A 



(9) 



The disjunction (V) of the first k variables enforces 
that at least one variable must be true, but the rest of 
{k{k — l)/2) 2-SAT type constraints ensure that only one 
of them is allowed to be true. 
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